The present invention relates to selecting optimal views for measuring radiation emitted by an object. More particularly, the present invention relates to selecting the optimal views using a scoring function based on information-theoretic measures. Of particular interest is view selection for medical imaging and/or in conjunction with medical instruments, such as guided minimally invasive surgical instruments.
Radionuclide imaging is one of the most important applications of radioactivity in medicine. Its purpose is to obtain a distribution image of a radioactively labeled substance, e.g., a radiopharmaceutical, within the body following administration thereof to a patient. Radioactive-emission imaging relies on the fact that in general, pathologies, such as malignant tumors, malfunctioning organs, and inflammations, display a level of activity different from that of healthy tissue. Thus, radiopharmaceuticals, which circulate in the blood stream, are picked up by the active pathologies to a different extent than by the surrounding healthy tissue; in consequence, the pathologies are operative as radioactive-emission sources and may be detected by radioactive-emission imaging. It will be appreciated that the pathology may appear as a concentrated source of high radiation, a hot region, as may be associated with a tumor, or as a region of low-level radiation, which is nonetheless above the background level, as may be associated with carcinoma.
A reversed situation is similarly possible. Dead tissue has practically no pick up of radiopharmaceuticals, and is thus operative as a cold region.
Thus radiopharmaceuticals may be used for identifying active pathologies as well as dead tissue.
In the discussion that follows, the term organ target is intended to include pathological features within organs. These pathological features may be expressed, by radioactive-emission imaging, as any one of the following:    i. hot regions, of a radioactive emission intensity higher than the background level;    ii. regions of low-level radioactive emission intensity, which is nonetheless above the background level; and    iii cold regions, of a radioactive emission intensity, lower than the background level.
Examples of radiopharmaceuticals include monoclonal antibodies or other agents, e.g., fibrinogen or fluorodeoxyglucose, tagged with a radioactive isotope, e.g., 99Mtechnetium, 67gallium, 201thallium, 111indium, 123iodine, 125iodine and 18fluorine, which may be administered orally or intravenously. The radiopharmaceuticals are designed to concentrate in the area of a tumor, and the uptake of such radiopharmaceuticals in the active part of a tumor, or other pathologies such as an inflammation, is higher and more rapid than in the tissue that neighbors the tumor. Thereafter, a radiation-emission-measuring-probe, which may be configured for extracorporeal or intracorporeal use, is employed for locating the position of the active area. Another application is the detection of blood clots with radiopharmaceuticals such as ACUTECT from Nycomed Amersham for the detection of newly formed thrombosis in veins, or clots in arteries of the heart or brain, in an emergency or operating room. Yet other applications include radioimaging of myocardial infarct using agents such as radioactive anti-myosin antibodies, radioimaging specific cell types using radioactively tagged molecules (also known as molecular imaging), etc.
The usual preferred emission for such applications is that of gamma rays, which emission is in the energy range of approximately 11-511 KeV. Beta radiation and positrons may also be detected.
Radioactive-emission imaging is performed with a radioactive-emission-measuring detector, such as a room temperature, solid-state CdZnTe (CZT) detector, which is among the more promising that is currently available. It may be configured as a single-pixel or a multi-pixel detector, and may be obtained, for example, from eV Products, a division of II-VI Corporation, Saxonburg Pa., 16056, or from IMARAD IMAGING SYSTEMS LTD., of Rehovot, ISRAEL, 76124, www.imarad.com, or from another source. Alternatively, another solid-state detector such as CdTe, HgI, Si, Ge, or the like, or a combination of a scintillation detector (such as NaI(Tl), LSO, GSO, CsI, CaF, or the like) and a photomultiplier, or another detector as known, may be used.
FIGS. 1a-1i schematically illustrate detecting units 102 and detecting blocks 101 of various geometries and constructions, and radioactive-emission-measuring probes associated with them.
FIG. 1a schematically illustrates a detecting unit 102, formed as a single-pixel detector 104, for example, a room-temperature solid-state CdZnTe (CZT) detector, having a diameter D and a thickness τd. Both the detector diameter D, or a diameter equivalent in the case of a non-circular detector, and the detector thickness τd affect the detecting efficiency. The detector diameter determines the surface area on which radioactive emission impinges; the greater the surface area, the greater the efficiency. The detector thickness affects the stopping power of the detector. High energy gamma rays may go through a thin detector, and the probability of their detection increases with detector thickness. By itself, a single-pixel detector cannot generate an image; rather, all counts are distributed over the surface area of the detector.
FIG. 1b schematically illustrates the detecting unit 102 with a collimator 108, formed as a single cell of a diameter D, a length L, and a septa thickness τ, attached to the detector 104. The collimator 108 may be, for example, of lead, tungsten or another material which substantially blocks gamma and beta rays.
The collimator's geometry, and specifically, the ratio of D/L, provides the detecting unit 102 with a collection angle δ analogous to a viewing angle of an optical camera. The collection angle δ limits the radioactive-emission detection to substantially only that radioactive emission, which impinges on the detector 104 after passing through a “corridor” of the collimator 108 (although in practice, some high-energy gamma rays may penetrate the collimator's walls).
FIG. 1c schematically illustrates a block 101 of the detecting units 102, with the collimator 108, formed as a multi-cell collimator, of a cell diameter D. The collection angle δ is defined for each of the detecting units 102 in the block, and each of the detecting units 102 forms a pixel in the block 101.
FIG. 1d schematically illustrates a radioactive-emission-measuring probe 100 which comprises several detecting units 102, of different geometries and different collection angles δ, within a housing 107.
FIGS. 1e-1i schematically illustrate the block 101, formed as a combination of a scintillation detector (such as NaI(Tl), LSO, GSO, CsI, CaF, or the like), a collimator grid, and photomultipliers.
As seen in FIG. 1e, the block 101, having proximal and distal ends 109 and 111, respectively, vis a vis an operator (not shown), is formed of the scintillation detector 104, of a single pixel, and the collimators 108, to create the detecting units 102. A plurality of photomultipliers 103 is associated with the single pixel scintillation detector 104, and with proper algorithms, as known, their output can provide a two dimensional image of the scintillations in the single pixel scintillation detector 104. In essence, this is an Anger camera, as known.
The distal view 111 of the collimator grid is seen in FIG. 1f. 
Two optional proximal views 109 of the photomultipliers 103 are seen in FIGS. 1g and 1h, as a square grid arrangement, and as an arrangement of tubes.
An Anger camera 117, of the block 101 in the housing 107 is seen in FIG. 1i. 
In each of the cases of FIGS. 1a-1i, the geometry of the collimator 108 determines the collection angle δ, wherein with no collimator, the collection angle δ, is essentially a solid angle of 4π steradians. Thus, the collimator's geometry affects both the detection efficiency and the image resolution, which are defined as follows:
i. The detection efficiency is the ratio of measured radiation to emitted radiation; and
ii. The image resolution is the capability of making distinguishable closely adjacent radioactive-emission organ targets, or the capability to accurately determine the size and shape of individual radioactive-emission organ targets.
Naturally, it is desired to optimize both the detection efficiency and the image resolution. Yet, they are inversely related to each other. The detection efficiency increases with increasing collimator's collection angle, and the image resolution decreases with increasing collimator's collection angle. For example, when the ratio of D/L is 1/2, the collection angle δ is substantially 2.5 steradians, so the cell views incident radiation within the confinement of about a 2.5-steradian sector. However, when the ratio of D/L is 1/12, the collection angle δ is substantially 0.31 steradians, so the cell views incident radiation within the confinement of about a 0.31-steradian sector.
Once the emission data is obtained, the data is processed to reconstruct the intensity distribution within the measured volume. The reconstruction process is generally complex, due to the large quantity of data which must be processed in order to obtain an accurate reconstruction. The following prior art statistical model may be used to perform reconstruction.
We assume an intensity distribution, I, defined over an input space U, where U comprises a set of basic elements (e.g., pixels in two dimensional spaces, voxels in three dimensional spaces), and I(u) is the intensity of a given basic element μ∈ U. A detecting unit positioned on a radiation-emission-measuring-probe takes a series of measurements y=(yt)t=1T from different positions and orientations around the volume U. The geometrical and physical properties of the detecting unit, together with its position and orientation in a given measurement t, determine the detection probability øt(u) of a photon emitted from location u. Thus the effective intensity of location u as viewed by the detecting unit during measurement t is øt(u)I(u).
The random count Xt(u) of photons that are emitted from location u and detected in measurement t is modeled by a Poisson process with mean øt(u)I(u). The total count of photons detected in measurement t is Yt=Σu∈U Xt(u), and the reconstruction problem is to reconstruct the intensities (I(u))u∈U from the measurements (yt)t=1T.
The 2-D Radon transform is a mathematical relationship which may be used to reconstruct the emission intensities of volume U when the set of measurements (yt)t=1T is unconstrained. The Radon transform is not statistical and does not take into account the Poissonian nature of the counts. In addition, it models the views as line projections. The Radon transform maps the spatial domain (x,y) to the Radon domain (p,φ). For a fixed projection angle, the Radon transform is simply a projection of the object. A technique known in the art as filtered back-projection (FBP) uses a back-projection operator and the inverse of the Radon transform to reconstruct the intensity distribution in volume U from measurements (yt)t=1T.
The basic, idealized problem solved by the FBP approach is to reconstruct an image from its Radon transform. The Radon transform, when properly defined, has a well-defined inverse. However, in order to invert the transform one needs measured data spanning 180°. In many medical imaging situations, the positioning of the detecting unit relative to the emitting object is constrained, so that complete measured data is not available. Reconstruction based on filtered back-projection is therefore of limited use for medical imaging. Maximum likelihood (ML) and Maximum A Posteriori (MAP) estimation methods, which address the statistical nature of the counts, have been found to provide better image reconstructions than FBP.
Limited-angle tomography is a reconstruction technique in the related art which reconstructs an image from projections acquired over a limited range of angular directions. The success of the reconstruction process depends upon the extent of the angular range acquired compared with the angular range of the missing projections. Any reconstruction from a limited range of projections potentially results in spatial distortions (artifacts) in the image. Limited angle techniques can be applied for both the Radon transform and the statistical models, but better results are generally achieved within the statistical framework. While it is known that the severity of the artifacts increases with the increasing angular range of the missing projections, limited-angle tomography does not provide information on which projections should be used in order to most effectively reconstruct the image.
Maximum likelihood (ML) estimation is a widely used method in the related art for reconstructing an image from a constrained set of measurements. A parameterization of the generative model described above is obtained by assigning an intensity I(u) to every voxel in U. The likelihood of the observed data y=(yt)t, given the set of parameters I={I(u):u∈U} is:
                              L          ⁡                      (                          y              ❘              I                        )                          =                              ln            ⁢                                                  ⁢                          P              ⁡                              (                                  y                  ❘                  I                                )                                              =                    ⁢                                    ln              ⁢                                                ∏                  t                                ⁢                                  P                  ⁡                                      (                                                                  y                        t                                            ❘                      I                                        )                                                                        =                                          ∑                t                            ⁢                              ln                ⁢                                                                  ⁢                                  P                  ⁡                                      (                                                                                            ∑                          u                                                ⁢                                                                              x                            t                                                    ⁡                                                      (                            u                            )                                                                                              ❘                      I                                        )                                                                                                              =                ⁢                              ∑            t                    ⁢                      ln            ⁢                                                  ⁢            Poisson            ⁢                                                  ⁢                          (                                                y                  t                                ❘                                                      ∑                    u                                    ⁢                                                                                    ϕ                        t                                            ⁡                                              (                        u                        )                                                              ⁢                                          I                      ⁡                                              (                        u                        )                                                                                                        )                                                              =                ⁢                              ∑            t                    ⁢                      {                                                                                                      -                                                                        ∑                          u                                                ⁢                                                                                                            ϕ                              t                                                        ⁡                                                          (                              u                              )                                                                                ⁢                          I                          ⁢                                                      (                            u                            )                                                                                                                +                                                                                                                                                                  y                        t                                            ⁢                      ln                      ⁢                                                                        ∑                          u                                                ⁢                                                                                                            ϕ                              t                                                        ⁡                                                          (                              u                              )                                                                                ⁢                                                      I                            ⁡                                                          (                              u                              )                                                                                                                                            -                                          ln                      ⁡                                              (                                                                              y                            t                                                    !                                                )                                                                                                                  }                              Note that the lower and upper bound of an indexing variable (such as voxels u and time index t) are omitted in the following description, when they are clear from the context.
There is currently no analytic way to solve Eqn. 1 for the maximum of the likelihood function. However, optimization methods that find local maxima of the likelihood are known. One such method is the Expectation-Maximization (EM) process.
Since the data generated by the model is only partially observable by our measurements, a basic ingredient of the Expectation-Maximization formalism is to define a set of random variables that completely define the data generated by the model. In the current case, since Yt=ΣuXt(u), the set of variables {Xu(t):U∈U; t=1, . . . , T} is such a set; the generated data is x=(xt)t, where xt=(xt(u))u, and the observed data y is completely determined by x. The main tool in the EM formalism is the complete data likelihood:
                                                                        ln                ⁢                                                                  ⁢                                  P                  ⁡                                      (                                          x                      ❘                      I                                        )                                                              =                            ⁢                              ln                ⁢                                                      ∏                    t                                    ⁢                                      P                    ⁡                                          (                                                                        x                          t                                                ❘                        I                                            )                                                                                                                                              =                            ⁢                                                ∑                  t                                ⁢                                  ln                  ⁢                                                            ∏                      u                                        ⁢                                          Poisson                      ⁢                                                                                          ⁢                                              (                                                                                                            x                              t                                                        ⁡                                                          (                              u                              )                                                                                ❘                                                                                                                    ϕ                                t                                                            ⁡                                                              (                                u                                )                                                                                      ⁢                                                          I                              ⁡                                                              (                                u                                )                                                                                                                                    )                                                                                                                                                                    =                            ⁢                                                ∑                  t                                ⁢                                                      ∑                    u                                    ⁢                                      {                                                                                            -                                                                                    ϕ                              t                                                        ⁡                                                          (                              u                              )                                                                                                      ⁢                                                  I                          ⁡                                                      (                            u                            )                                                                                              +                                                                                                    x                            t                                                    ⁡                                                      (                            u                            )                                                                          ⁢                                                  ln                          ⁡                                                      (                                                                                                                            ϕ                                  t                                                                ⁡                                                                  (                                  u                                  )                                                                                            ⁢                                                              I                                ⁡                                                                  (                                  u                                  )                                                                                                                      )                                                                                              +                                              ln                        ⁡                                                  (                                                                                                                    x                                t                                                            ⁡                                                              (                                u                                )                                                                                      !                                                    )                                                                                      }                                                                                                          (        2        )            
Since the likelihood depends on the complete data, which is only partially observable, we take its expectation with respect to the space of the unobserved data, given the current set of hypothesized parameters (i.e. the current estimator). The result is a function Q(I|I′) which assigns likelihood to sets I of model parameters, given the current set I′, and given the observed data y:
                                                                        Q                ⁡                                  (                                      I                    |                                          I                      ′                                                        )                                            =                            ⁢                              E                ⁡                                  [                                                            ln                      ⁢                                                                                          ⁢                                              P                        (                                                  x                          ⁢                                                                                  I                            )                                                                                                                      ⁢                      y                                        ;                                          I                      ′                                                        ]                                                                                                        =                            ⁢                                                ∑                  t                                ⁢                                                      ∑                    u                                    ⁢                                      {                                                                                            -                                                                                    ϕ                              t                                                        ⁡                                                          (                              u                              )                                                                                                      ⁢                                                  I                          ⁡                                                      (                            u                            )                                                                                              +                                              E                        ⁡                                                  [                                                                                                                                                      x                                  t                                                                ⁡                                                                  (                                  u                                  )                                                                                            |                                                              y                                t                                                                                      ;                                                          I                              ′                                                                                ]                                                                                                                                                                                                          ⁢                                                ln                  ⁡                                      (                                                                                            ϕ                          t                                                ⁡                                                  (                          u                          )                                                                    ⁢                                              I                        ⁡                                                  (                          u                          )                                                                                      )                                                  +                C                            }                                                          (        3        )            where C is a term which is independent of the intensities I. The function Q(I|I′) is maximized by the following new estimates:
                                          I            ⁡                          (              u              )                                =                                    1                                                ∑                  t                                ⁢                                                      ϕ                    t                                    ⁡                                      (                    u                    )                                                                        ⁢                                          ∑                t                            ⁢                              E                ⁡                                  [                                                                                                              x                          t                                                ⁡                                                  (                          u                          )                                                                    |                                              y                        t                                                              ;                                          I                      ′                                                        ]                                                                    ;                  ∀                      u            ∈                          U              .                                                          (        4        )            
The expectation in Equation 4 is obtained as follows:
                                                                        P                                                      X                    t                                    ⁡                                      (                    u                    )                                                              (                                                                                          x                      t                                        ⁡                                          (                      u                      )                                                        ⁢                                                                                                        y                        t                                            ;                                              I                        ′                                                              )                                                  =                                                                                                    P                                                  Y                          t                                                                    ⁡                                              (                                                                                                            y                              t                                                        |                                                                                          x                                t                                                            ⁡                                                              (                                u                                )                                                                                                              ;                                                      I                            ′                                                                          )                                                              ⁢                                                                  P                                                                              X                            t                                                    ⁡                                                      (                            u                            )                                                                                              ⁡                                              (                                                                                                            x                              t                                                        ⁡                                                          (                              u                              )                                                                                |                                                      I                            ′                                                                          )                                                                                                                        P                                              Y                        t                                                              ⁡                                          (                                                                        y                          t                                                |                                                  I                          ′                                                                    )                                                                                                                                              =                                                                                                                  Poisson                        (                                                                                                            y                              t                                                        -                                                                                          x                                t                                                            ⁡                                                              (                                u                                )                                                                                                              |                                                                                    ∑                                                              v                                ≠                                u                                                                                      ⁢                                                                                                                            ϕ                                  t                                                                ⁡                                                                  (                                  v                                  )                                                                                            ⁢                                                                                                I                                  ′                                                                ⁡                                                                  (                                  v                                  )                                                                                                                                                                    )                                                                                                                                                Poisson                        ⁡                                                  (                                                                                                                    x                                t                                                            ⁡                                                              (                                u                                )                                                                                      |                                                                                                                            ϕ                                  t                                                                ⁡                                                                  (                                  u                                  )                                                                                            ⁢                                                                                                I                                  ′                                                                ⁡                                                                  (                                  u                                  )                                                                                                                                              )                                                                                                                                      Poisson                  (                                                            y                      t                                        |                                                                  ∑                        v                                            ⁢                                                                                                    ϕ                            t                                                    ⁡                                                      (                            v                            )                                                                          ⁢                                                  I                          ⁡                                                      (                            v                            )                                                                                                                                )                                                                                                        =                              Binomial                (                                                                                                    x                        t                                            ⁡                                              (                        u                        )                                                              |                                                                                                                        ϕ                            t                                                    ⁡                                                      (                            u                            )                                                                          ⁢                                                                              I                            ′                                                    ⁡                                                      (                            u                            )                                                                                                                                                ∑                          v                                                ⁢                                                                                                            ϕ                              t                                                        ⁡                                                          (                              v                              )                                                                                ⁢                                                                                    I                              ′                                                        ⁡                                                          (                              v                              )                                                                                                                                                            ;                                      y                    t                                                  )                                                                        (        5        )            
It follows that
            E      [                                                  x              t                        ⁡                          (              u              )                                |                      y            t                          ;                  I          ′                    ]        =                  y        t            ⁢                                    ϕ            t                    ⁢                      (            u            )                    ⁢                                    I              ′                        ⁡                          (              u              )                                                            ∑            v                    ⁢                                                    ϕ                t                            ⁡                              (                v                )                                      ⁢                                          I                ′                            ⁡                              (                v                )                                                          ,and hence the EM iteration is:
                              I          ⁡                      (            u            )                          =                              1                                          ∑                t                            ⁢                                                ϕ                  t                                ⁡                                  (                  u                  )                                                              ⁢                                    ∑              t                        ⁢                                          y                t                            ⁢                                                                    ϕ                    t                                    ⁢                                      (                    u                    )                                    ⁢                                                            I                      ′                                        ⁡                                          (                      u                      )                                                                                                            ∑                    v                                    ⁢                                                                                    ϕ                        t                                            ⁡                                              (                        v                        )                                                              ⁢                                                                  I                        ′                                            ⁡                                              (                        v                        )                                                                                                                                                    (        6        )            
It is provable that each EM iteration improves the likelihood. Thus, given a random starting estimator, the EM algorithm iterates the above improvement step until it converges to a local maximum of the likelihood. Several random starts increase the chance of finding a globally good estimator.
It is usually desired to maximize the expected posterior probability (given a proper prior) rather than the expected likelihood. In that case we assume a prior probability on the intensities P(I)=πuP(I(u)). A proper conjugate prior for the Poisson distribution is the Gamma distribution:
                              P          ⁡                      (                          I              ⁡                              (                u                )                                      )                          =                              Gamma            (                                                            I                  ⁡                                      (                    u                    )                                                  |                                  α                  u                                            ;                              β                u                                      )                    ⁢                                    β              u                                                α                  u                                +                1                                                    Γ              ⁡                              (                                                      α                    u                                    +                  1                                )                                              ⁢                                                    I                ⁡                                  (                  u                  )                                            ⁢                                                                                  α              u                                ⁢                      ⅇ                                          -                                  β                  u                                            ⁢                              I                ⁡                                  (                  u                  )                                                                                        (        7        )            
Now the maximization is done on Q(I|I′)=E[lnP(x|I)p(I)|y;I′]. Plugging the Gamma prior into Q, and solving for I(u), we get the following EM iteration for the maximum posterior estimation:
                                                                        I                ⁡                                  (                  u                  )                                            =                            ⁢                                                                    α                    u                                    +                                                            ∑                      t                                        ⁢                                          E                      ⁡                                              [                                                                                                                                            x                                t                                                            ⁡                                                              (                                u                                )                                                                                      |                                                          y                              t                                                                                ;                                                      I                            ′                                                                          ]                                                                                                                                  β                    u                                    +                                      ∑                                                                  ϕ                        t                                            ⁡                                              (                        u                        )                                                                                                                                                                    =                            ⁢                                                1                                                            β                      u                                        +                                                                  ∑                        t                                            ⁢                                                                        ϕ                          t                                                ⁡                                                  (                          u                          )                                                                                                                    [                                                      α                    u                                    +                                                            ∑                      t                                        ⁢                                                                  y                        t                                            ⁢                                                                                                    ϕ                            t                                                    ⁢                                                      (                            u                            )                                                    ⁢                                                                                    I                              ′                                                        ⁡                                                          (                              u                              )                                                                                                                                                            ∑                            v                                                    ⁢                                                                                                                    ϕ                                t                                                            ⁡                                                              (                                u                                )                                                                                      ⁢                                                                                          I                                ′                                                            ⁡                                                              (                                v                                )                                                                                                                                                                                                        ]                                                                                                            (              8              )                                                                          (              9              )                                          
The EM update step can be formulated in matrix notation as follows. Let Φ be the matrix of the projections [φt(u)]t,u, and let I,I′,y, α and β be represented as column vectors. Equation 8 can be written in vector and matrix notations as:
                    I        =                              α            +                                          I                ′                            ·                              (                                                      Φ                    T                                    ⁢                                      y                                          Φ                      ⁢                                                                                          ⁢                                              I                        ′                                                                                            )                                                          β            +                                          Φ                T                            ⁢              1                                                          (        10        )            where the explicit multiplication and division denote element-wise operations, and where 1 is a vector (of the appropriate length) consisting solely of 1′s.
Limited computational resources (i.e., when the entire projection matrix Φ cannot be kept in memory) may require breaking the update computation according to a partition of Φ into a set of sub-matrices (Φi). In that case the intensities can be updated gradually (using only one sub-matrix at each step) according to the following computation:
                    I        =                              α            +                                          I                ′                            ·                                                ∑                  i                                ⁢                                                      ϕ                    i                    T                                    ⁢                                                            y                      t                                                                                      Φ                        t                                            ⁢                                              I                        ′                                                                                                                                      β            +                                          ∑                i                            ⁢                                                Φ                  i                  T                                ⁢                1                                                                        (        11        )            where yi is the vector of observations that are obtained using the views of Φi.
In order to achieve a reconstructed image which is adequate for medical diagnostic and treatment purposes, a high-resolution image of the tested object must be obtained. When high-resolution detecting units are used, their efficiency is relatively low, and the detecting units must remain at each position for a relatively long time in order to achieve a high probability of detection. Since during medical testing, measurements are generally performed at many locations as the detecting unit is moved relative to the observed organ, the testing procedure generally requires a long time and is physically and emotionally difficult for the patient. Additionally, reconstruction is based upon a large quantity of data, and is a lengthy and computationally complex process.